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Abstract. The global population structure and dispersal patterns of Pantala flavescens (Fa- 
bricius, 1798) are evaluated using a geographically extensive mitochondrial DNA dataset, a 
more limited samples of nuclear markers, wing isotopic (52H) data and a literature review. 
No spatial or temporal haplotype structure was recovered between the samples. Isotope 
data suggest that most samples were immigrants at the collection locations. A literature 
review of migration events for the species confirms regular inter-and intra-continental mi- 
grations occur (the majority reported from Asia, Africa and Australasia), with individuals 
and swarms dispersing thousands of kilometers over land and oceans. Migrations coincide 
with prevailing winds and seasonal rains, which points to a mechanism we name the “pan- 
tropical Pantala conveyor belt”, suggesting widespread gene flow is possible for an aquatic 
insect with excellent flying ability linked to rapid larval development. 
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Introduction 


The Globe Skimmer or Wandering Glider dragonfly, Pantala flavescens (Fabricius, 
1798) (Odonata: Anisoptera: Libellulidae), is a well-known long-distance migrant. 
It has a circum-global distribution, and is the most widely distributed of all dragon- 
fly species (Fraser, 1936; Russell et al., 1998). Pantala flavescens occurs most often 
in warmer climates, especially within the tropics, but also migrates seasonally into 
temperate zones (Borisov & Malikova, 2019; May, 2013). Its ability to disperse over 
great distances is demonstrated by its occurrence on offshore islands, for example 
on many Pacific Ocean islands (Rowe, 2004; Schmidt, 1938) including New Zealand 
(Corbet, 1979; Lieftinck, 1975; Rowe, 1980). Pantala flavescens has even reached 
the remote Sub-Antarctic Amsterdam Island (Devaud & Lebouvier, 2019) and is 
established on Rapa Nui (formerly Easter Island) (Dumont & Verschuren, 1991; 
Moore, 1993; Samways & Osborn, 1998). Its longest, regular, well-documented 
migratory route is across the western Indian Ocean (Anderson, 2009; Hobson et 
al., 2012a), where it has been suggested to follow seasonal (monsoonal) rains, 
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which allow females to take advantage of temporary 
pools of water for oviposition (Anderson, 2009; Corbet, 
1999; Suhling et al., 2015). 

The ability of P. flavescens to migrate such long dis- 
tances is facilitated by a unique combination of morpho- 
logical characteristics. These include distinct wing mor- 
phology (Alvial et al., 2019; Li et al., 2014; Moore, 1993; 
Sacchi & Hardersen, 2013; Suarez-Tovar & Sarmiento, 
2016; Outomuro & Johansson, 2019; Zhao, 2012), of 
which the most obvious adaptation is their enlarged 
hind wing bases (a feature associated with gliding, 
Corbet, 1999), as well as unique thoracic musculature 
(Bäumler et al., 2018). Physiological features that allow 
individuals to migrate long distances likely involve their 
ability to make efficient use of lipids as fuel rather than 
carbohydrates (e.g., Kallapur & George, 1973); it has 
been suggested that dragonfly fat reserves can sustain 
eight hours of flapping flight. In addition, P. flavescens 
may be able to feed on insects while migrating thereby 
replenishing their fat stores en route (Anderson, 2009). 

Given its wide distribution and extraordinary migra- 
tory abilities, there is interest in the population struc- 
ture of P. flavescens. Studies in single countries have 
revealed high haplotype diversity and an expected lack 
of population structure (Cao et al., 2015; Hayashi & 
Arai, 2004; Low, 2017). Over a much wider geographi- 
cal range, albeit with a small sample size, Troast et al. 
(2016) found that P. flavescens specimens from South 
America, North America, and Asia shared mitochondrial 
DNA haplotypes, suggesting high rates of gene flow and 
a single global panmictic population. Similarly, Alvial et 
al. (2017) found a lack of population structure between 
samples from across South and Central America using 
both mitochondrial and nuclear DNA markers. In con- 
trast, a subsequent study by Alvial et al. (2019), which 
compared samples not only from Central and South 
America but also from islands in the Pacific and Indian 
Oceans, did find some suggestion of population struc- 
ture between ‘continental’ and ‘insular’ samples using 
microsatellite loci, which also suggested that individuals 
from Rapa Nui formed a genetically distinct population. 
In addition, Pfeiler & Markow (2017) reanalysed pub- 
lished data (Low, 2017; Troast et al. 2016) and suggest- 
ed that ‘concluding global panmixia in P. flavescens may 
be premature.’ Thus, the results from previous studies 
are ambiguous, with some suggesting panmixia while 
others suggest some degree of population structure. 

Global panmixia, or random mating and uninterrupt- 
ed gene flow across an entire species, appears to be un- 
common in widespread species. Barriers to genetic ex- 
change, which may be geographical, temporal or behav- 
ioural in nature, tend to lead to the evolution of geneti- 
cally distinct regional populations. However, in some 
species dispersal or migration can lead to widespread 
gene flow despite the presence of obstacles (Neeth- 
ling et al., 2008; Oomen et al. 2011). A commonly cited 
example is that of the European Eel (Anguilla anguilla 
(Linnaeus, 1758)), in which individuals from separate 
river basins travel to the Sargasso Sea for reproduction, 
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leading to panmixia (Als et al., 2011; Palm et al., 2009; 
but see Wirth et al., 2001 for counter argument). Wide- 
spread genetic mixing has also been recorded in several 
other groups of organisms, including insects, birds, bats 
and plants (Als et al., 2001; Naro-Maciel, 2011; Neeth- 
ling et al., 2008; Oomen et al., 2011; Ridgway et al., 
2001; Peel et al., 2013; White et al., 2011). Neverthe- 
less, most reports of panmixia might best be described 
as ‘local panmixia’, where random mating occurs across 
a population in a particular region rather than across 
the entire species (e.g. Peel et al., 2013; White et al., 
2011). Studies on the population structure of long- 
range migrant insects are still rare, with the exception 
of locusts (Locusta migratoria, (Linnaeus, 1758) and 
Schistocerca gregaria Forsskal, 1775), monarch butter- 
flies (Danaus plexippus (Linnaeus, 1758)), painted lady 
butterflies (Vanessa cardui (Linnaeus, 1758)), and the 
green darner dragonfly, Anax junius (Drury, 1773) (e.g., 
Hallworth et al., 2018). 

Here, we evaluate the global population structure of 
P. flavescens, using genetic data to test the hypothesis of 
mixing across its global range. We have compiled genetic 
sequences available through the NCBI GenBank online 
database, supplementing this with a limited number of 
new samples, to form the largest such dataset yet as- 
sembled, including samples from 29 countries across six 
continents. We estimated various population structure 
measures, including testing for levels of gene flow. Fur- 
ther, we measured stable hydrogen isotope (57H) values 
in a small sample of P. flavescens wings to establish what 
proportion of these specimens originated in regions dif- 
ferent from those where they were collected (Hayashi 
& Arai, 2004; Hobson et al., 2012a, 2012b; 2021). This 
allowed us to test for the degree of migratory behaviour 
in P. flavescens, on the assumption that the presence of 
a large number of immigrant individuals across a global 
range would confirm widespread migration. Finally, we 
compiled published observations of migratory behav- 
iour which served to identify possible connections and 
pathways between and within the continents. 


Materials and methods 
Study organism 


Pantala flavescens (Fabricius, 1798) is a member of the 
dragonfly family Libellulidae (Odonata: Anisoptera). 
The vernacular English names “globe skimmer” and 
“wandering glider” reflect this species’ long-distance 
flying ability and wide distribution. Adults are red, yel- 
low, or brown in colour, commonly observed “hawking” 
over open areas. Egg-laying often occurs in tandem, 
and females produce their eggs continuously (Koch et 
al., 2011; Sharma, 2017; Ware et al., 2012). Eggs are 
distributed widely, even being oviposited among differ- 
ent ponds during a single oviposition event, which is in- 
terpreted as risk-spreading to avoid loss of all offspring 
should temporary habitats dry out before development 
is completed (Schenk et al., 2004). Larvae mostly inhab- 
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it temporary or ephemeral lentic waters including rain- 
pools forming after heavy rainfalls. Offspring develop 
very rapidly, taking 4-6 weeks from egg hatch to adult 
emergence (e.g., Hawking & Ingram, 1994; Kumar, 
1984; Suhling et al., 2004). 


Sampling 


We collected 23 specimens of P. flavescens from six 
countries during 2004 to 2015 (see Supplementary ma- 
terial, Table $1). In addition we downloaded genetic se- 
quences of P. flavescens from NCBI GenBank (N = > 600 
for COI; complete list in Supplementary material, Table 
S1). We included COI sequences from North America 
(N = 16; Canada, USA), Central and South America (N = 
25; Brazil, Chile, Costa Rica, Guyana, Peru), Asia (N 
190; Cambodia, China, India, Indonesia, Japan, Korea, 
Malaysia, Maldives, Pakistan, Philippines, Singapore, 
Vietnam), Africa (N = 4; Guinea Bissau, Liberia, Sen- 
egal; unfortunately, samples from South Africa and 
Namibia failed during the sequencing step), Australia 
and the Pacific (N = 13; Rapa Nui [a territory of Chile], 
Fiji, French Polynesia, Hawaii, Mariana Islands, Tonga 
Islands). Additionally, we downloaded and ran analy- 
ses on P. flavescens nuclear genetic data available from 
NCBI (Table S1). 


Genetic dataset partitions 


We ran our analyses with the geographic dataset de- 
scribed above for the COI. Additionally, we ran the COI 
haplotype analyses with either (1) fewer Japanese sam- 
ples or (2) excluding samples from Alvial et al. (2017, 
2019). 

1. Inference of population structure can be sensitive to 
sampling bias, (see Puechmaille, 2016). More than 
half of our samples came from Japan, while other re- 
gions had < 10 individuals. Thus, for additional hap- 
lotype analyses we limited the number of Japanese 
samples to 77 (“reduced” data set). We chose these 
samples to sample P. flavescens evenly from across 
the Japanese islands for which we had data. This 
gave us 248 P. flavescens sequences in total. 

2. For NCBI samples from Alvial et al.’s (2017, 2019 
studies (N = 40)) it was impossible to distinguish 
which samples were which. Hence, these were 
coded as “Alvial” (i.e. grouping together all of her 
samples from Costa Rica, Chile, Peru, Tonga, Rapa 
Nui, and the Maldives). Alvial et al.’s samples were 
included in testing our Null hypothesis that there is 
no global geographic structure. These samples had 
to be excluded from the rest of our geographic hy- 
pothesis testing scenarios as we could not infer prov- 
enance (i.e., if the samples were from the Northern 
or Southern, Eastern, or Western Hemisphere, see 
Inference of Population Structure section for more 
details on scenarios). 
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We additionally reconstructed haplotype networks 
for P. flavescens sequences available from NCBI Gen- 
Bank for mitochondrial 16S and Cytochrome B, as well 
as nuclear 18S and Histone 3. Apart from 16S, these 
sequence datasets were from single countries, and so 
although they do provide insights into genetic diversity 
they shed no direct light on global population struc- 
ture. 


DNA Extraction/PCR Amplification 


DNA was extracted from tissue of P. flavescens samples 
which had been stored in ethanol or acetone-dried. 
Extraction was completed using Qiagen DNeasy Blood 
and Tissue Kits (Maryland, USA). We sequenced the 
mitochondrial Cytochrome Oxidase | (COI; ~ 658 base 
pairs), a gene which was chosen because of a large 
number of existing sequences available on NCBI, allow- 
ing for this large, global dataset. Polymerase Chain Re- 
action (PCR) was done using the primer set of Simon 
et al. (1994): 1709f (TAATTGGAGGATTTGGAAATTG, T „= 
56.2°C) and 2191r (CCYGGTARAATTARAATRTARACTTC, 
T, =56.7°C). PCR reaction total volume was 25 uL, com- 
prising 12.5 uL Taq 2X Master Mix, 5-7 uL DNA, 3.5- 
5.5 uL RNA-ase free water, and 1-2 uL of the forward 
and reverse primers. 

PCR was done using an Eppendorf Thermocycler us- 
ing the following parameters: initial heat step at 94.0°C 
held for one minute, 20 cycles of denaturation (30 s, 
94.0°C), annealing (45 s, 48.0°C), and extension (45 s, 
72.0°C), an additional 20 cycles of denaturation (45 s, 
94.0°C), annealing (45 s, 50.0°C), and extension (45 s, 
72.0°C), followed by a final extension at 72.0°C held for 
five minutes. PCR products were visualized via gel elec- 
trophoresis. We used commercially available Macrogen 
services (New York City, New York, USA) for sequencing 
and purification. 


Sequences 


Sequences were edited using Sequencher® (2016). 
These sequences, along with those from NCBI Gen- 
Bank, were aligned using Clustal-X Version 2.1 (Larkin, 
2007) and Mesquite Version 2.75 (Maddison & Maddi- 
son, 2017). For the COI, twenty-three de novo sequenc- 
es from Africa, Australia, the Maldives, and the USA 
were sequenced and resulting sequences were depos- 
ited in NCBI GenBank (Accession numbers OM945847- 
OM945869). In total we had 25 samples from Central 
and South America, 16 from North America, 13 from 
Australasia, 190 from Asia including the Maldives, and 
four from Africa for the COI. For the 16S the 16 sequenc- 
es were from Africa, India, Japan, Korea, Malaysia, Paki- 
stan, and the United States of America, for the 18S the 
four sequences were from Japan, for the Cytochrome B 
the 77 sequences were from China, and for the Histone 
3 the four sequences were from Japan. 
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Haplotype analyses 


To assess variation among individuals of P. flavescens, 
we reconstructed minimum spanning networks for 
each gene fragment using PopArt version 1.7 (Leigh & 
Bryant, 2015) with the minimum spanning networks al- 
gorithm (Bandelt et al., 1999). 


Inference of population structure 


We also assessed population structure summary sta- 
tistics using the COI genetic data. All statistics were 
estimated from the “reduced dataset”, and required 
assigning individuals to a fixed number of putative geo- 
graphically based hypotheses (populations), which we 
did in several ways. First, we estimated the number of 
genetic clusters/populations using K-means clustering 
analysis in adegenet R package 2.1.3 (Jombart, 2008; 
Jombart & Ahmed 2011). We selected the optimal num- 
ber of clusters using the Bayesian Information Criterion 
(BIC); choosing the solution with the best BIC value. Dis- 
criminant Analysis of Principal Components (DAPC), as 
implemented in adegenet R package, was then used to 
assess genetic variation between identified population 
clusters. We also determined membership probability 
of each sample for each of the predicted populations; 
to do so we calculated the degree of genetic differentia- 
tion (F,,) among predicted clusters using the hierfstat R 
package version 0.5-7 (Goudet, 2005), and we ran an 
AMOVA in the poppr package version 2.8.6 (Kamvar et 
al., 2014, 2015) to test whether the observed variation 
was higher between or within each cluster. 

We estimated population genetics statistics assuming 
different scenarios considering possible biogeographi- 
cally separated populations: (a) assuming the popula- 
tions were defined based on the results from K-means 
clustering, (b) assuming each of the continents is a sepa- 
rate population, (c) assuming Northern/Southern Hemi- 
spheres are separate populations, (d) assuming Eastern/ 
Western Hemispheres are separate populations, and (e) 
a null hypothesis of one panmictic population. Insects 
generally do not follow geo-political boundaries, but we 
also ran an additional scenario (f) assuming countries 
are separate populations. Population genetics statistics 
under the various scenarios were estimated in the soft- 
ware package R as described above. We calculated the 
genetic structure index, Fa and other statistics (haplo- 
type diversity, Hd; nucleotide differences, K; number of 
polymorphic sites, S; Tajima’s D, Fu and Li’s statistics, 
and nucleotide diversity, m) under each of these popula- 
tion scenarios as described above. For the single pan- 
mictic population global scenario we could not estimate 
F as there was only one population. 


Effective migration rates 


For the COI, we used the program Migrate-N version 
3.7.2 (Beerli et al., 2019) to estimate effective migra- 
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tion rates (4N M) and © = 4N „u under scenarios (c) and 
(d). We let the migration matrix stay at default settings. 
We ran 10 replicate maximum likelihood Migrate-N 
runs of 20,000 steps per chain each with a Brownian 
motion model, and one long chain. We sampled every 
100 steps, with a burn-in of 10,000 steps and a bound- 
ed-adaptive heating scheme with four temperatures 
(3 cold, 1 hot) of 1.00, 1.50, 3.00, and 10,000. To test 
the effect of temperature on our results, we ran an 
analysis with two different “hottest temperatures” and 
the results did not differ. We evaluated the best model 
using Log Bayes Factor (LBF) and highest model prob- 
ability values in Migrate-N. Effective Sample Sizes (ESS) 
were evaluated in TRACER version 1.7.1 (Rambaut et 
al., 2018). 


Stable Isotope Analysis 


We excised wings from 29 individuals representing 
each of the continents of North America, Central and 
South America, Africa, Australia, and Asia plus the 
Pacific; samples were chosen from both islands and 
mainland areas (Table S2). All samples had been dried 
in acetone, which has been shown to have no impact 
on isotopic results (Hobson et al., 2012a; Hobson & 
Wassenaar, 2019). We submitted wings to the Cornell 
Isotope Laboratory (COIL) facility, where surface oils 
were removed using a 2:1 chloroform:methanol solu- 
tion and air dried in a fume hood. Subsamples of wings 
were weighed into silver capsules, crushed, and loaded 
into a zero-blank carousel attached to TC/EA (Thermo, 
Bremen, Germany) where samples were combusted at 
1350°C pyrolytically. Resultant H, gas was analyzed un- 
der continuous-flow isotope-ratio mass spectrometry 
(CFIRMS) on a Thermo Delta-V mass spectrometer us- 
ing a helium carrier gas. Samples were calibrated us- 
ing within-run measurements of keratin standards CBS 
(-197%o) and KHS (-54.1%o) according to the compara- 
tive equilibration method of (Wassenaar & Hobson, 
2003). An in-house (Spectrum) commercial standard 
was used to correct for instrument drift. Stable-hydro- 
gen isotope values are reported here in standard delta 
notation in parts per thousand (%o) deviation from the 
Vienna Mean Ocean Water (VSMOW) standard. Based 
on within-run replicate measurements of standards, we 
estimated measurement precision to be ~ 2%o. 
Assigning P. flavescens to geographical origins using 
the stable-hydrogen isotope approach is problematic 
because this species is known to use ephemeral pools 
that may be created from highly seasonal (often mon- 
soonal) rainfall events. Hobson et al. (2012a) modeled 
the origins of individuals captured on the Maldives en 
route to Africa using an amount-weighted mean annual 
precipitation model that provided the first isotopic evi- 
dence for long-distance migration in this species and 
pointed to origins in northern India or possibly even 
further north and east. That approach was based on a 
calibration algorithm derived from known-origin drag- 
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onflies in North America (Hobson et al., 2012b: Wing 
67H = 0.91(MAD) -42.5, r? = 0.75, where MAD is the 
amount-weighted mean annual precipitation for a giv- 
en site). Based on the findings of Lopez-Calderon et al. 
(2019), we used this relationship to evaluate whether 
a given dragonfly likely originated from sites of capture 
using the normal probability density function and with 
an assumed standard deviation of the residuals of the 
calibration algorithm to be a highly conservative 19%o 
(Hobson et al., 2012b). We considered three odds ratios 
(2:1, 3:1, and 4:1) for each individual assignment as be- 
ing a “local” or an “immigrant”. 


Results 
Haplotypes 


Our COI dataset consisted of an alignment with 686 
characters. Haplotype minimum spanning analyses on 
all datasets suggest that most P. flavescens samples 
share a single common haplotype for COI (Figure 1, net- 
work from reduced dataset); we found one main hap- 
lotype with 24—43 additional haplotypes varying by < 3 
nucleotides (43 = full dataset and 24 = reduced dataset). 
Regardless of what level of population assignment was 
used (continent, Western/Eastern Hemisphere, North- 
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ern/Southern Hemisphere), we recovered the same 
haplotype structure. For the other gene fragments stud- 
ied, one or two main haplotypes were recovered (Fig- 
ure S3), with the exception of Cytochrome B, for which 
several haplotypes were suggested among the samples 
(all of which were from China). In the case of the mi- 
tochondrial marker 16S, we located 16 samples from 
three continents on NCBI GenBank, all of which shared 
a single haplotype (Figure S3). In the case of the nuclear 
marker H3, only four sequences were available on NCBI 
Genbank, three of which shared the same haplotype; 
the fourth sample revealed a very different haplotype, 
and likely represents a different species, mislabeled as 
P. flavescens, or contamination. 


Population Genetics 


The clustering method in the adegenet package esti- 
mated K = 1 to three populations of P. flavescens as most 
likely for COI. We used the “elbow” method to choose 
K = 2 and K = 3 as the population cluster values (sce- 
nario (a); Jombart & Collins, http://adegenet.r-forge.r- 
project.org/files/tutorial-genomics.pdf). However, as 
inferred clusters higher than K 2 2 lacked any geograph- 
ic structure, we consider these values as to be difficult 
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to interpret: they do not seem to share any clear pat- 
tern that unites the individuals in these clusters (i.e., 
these clusters do not reflect timing of sampling, geo- 
graphic region of sampling, etc.). In addition, the BIC 
scores did not vary greatly among the K = 1, K = 2, and 
K = 3 (BIC score range 300—~450) and DAPC scatterplots 
showed large overlap across clusters (Figure 2). We fur- 
ther tested whether differences between K = 2 and K = 
3 could be explained by temporal bias in sampling, but 
as each putative K population contained samples from 
across sampling periods (i.e., early 2000s and early/mid 
2010s), we reject this explanation for inferred patterns. 
Lastly, under the K = 2 and K = 3 scenarios we found low 
genetic structure and F_ values between the popula- 
tions, consistent with high rates of gene flow (most F, 
values were < 0.5). 

F.. values between different biogeographical sce- 
narios varied, ranging from below 0.05 in scenario (b) 
(between Asia and Africa, Asia and Australasia, North 
America and Asia) indicating little population differen- 
tiation, to 0.55 between Africa and Australasia (Table 1) 
suggesting greater population differentiation. However, 
with only four samples from Africa, all of which were 
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from West Africa, this result should be treated with 
caution. For scenario (c), F, was 0.007; for scenario (d) 
F „ was 0.004. In short, for the geographical scenarios 
(b)-(d) our dataset suggests high gene flow and no ge- 
netic structure between putative populations. 


Effective migration rates 


Migrate-N found widespread gene flow in COI. Mi- 
grate-N analytical results suggested that individuals 
migrating (in either direction) from Western Hemi- 
sphere to Eastern Hemisphere was the best fit with 
In(Prob(D| Model)) = -4158.48 compared to the model 
for individuals migrating (in either direction) from the 
Northern/Southern Hemisphere (Northern/Southern 
Hemisphere In(Prob(D|Model)) = -6480.57). Within 
the Eastern/Western Hemisphere analysis, movement 
in both directions had normal distributions of migra- 
tion rates, M, although relative variance around the 
means varied. Mean migration rate for the Western 
Hemisphere-Eastern Hemisphere direction was 478.2, 
and the mean migration rate in the Eastern Hemisphere 
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Table 1. Pairwise F values under various population scenarios: (a) the K-means clustering analysis populations estimated in 
adegenet R package, (b) assuming continents are separate populations. 

















Scenario (a) K=2 K=3 

Population F, 2 1 3 2 1 

3 - = 0.0 0.217 0.158 
2 0 0.168 0.217 0 0.277 
1 0.168 0 0.158 0.277 0 
Within population variation 58.983 56.462 

Among population variation 41.017 43.538 

Scenario (b) 

Population F.. (b) Africa Asia Australasia N America S America 
Africa (N=4) 0 0.029 0.551 0.197 0.223 
Asia (N=190) 0.029 0 0.002 0.031 0.062 
Australasia (N=13) 0.551 0.002 0 0.158 0.111 
North America (N=16) 0.197 0.031 0.158 0 0.131 
South America (N=25) 0.223 0.062 0.111 0.131 0 
Within population variation 95.889 

Among population variation 4.111 





Table 2. Population genetics statistics based on different population scenarios, with K: nucleotide differences, m: nucleotide 


diversity, S: # of polymorphic sites. p < 0.05 in bold. 





Population Tajima’s D K T S Hd variance FuandLisD Fu and Li’s F Fst 
Western Hemisphere -1.37923 0.300 0.001 3 0.2360.009 -1.50777 -1.70431 0.004 
(N=41) p>0.10 p>0.10 p>0.10 
Eastern Hemisphere -2.18914 0.943 0.006 24 0.5730.002 -3.24440 -3.38651 
(N=207) p<0.01 p<0.05 p<0.02 
Northern Hemisphere -2.16518 0.494 0.003 20 0.5810.002 -3.69480 -3.69542 0.007 
(N=228) p<0.01 p<0.02 p<0.02 
Southern Hemisphere -1.79423 0.251 0.002 3 0.3280.009 -1.64715 -1.97431 
(N=20) 0.10>p>0.05 p>0.10 0.10>p>0.05 
All data -2.35627 0.452 0.003 22 0.3480.002 -4.45510 -4.34133 
p<0.01 p<0.02 p<0.02 





to Western Hemisphere direction was 40.0; the mean 
migration rate in the North-South direction was 311.1, 
and the mean migration rate in the South-North direc- 
tion was 19.5. In short, the results of this analysis sug- 
gest more migration going from west to east and north 
to south than in the opposite directions. 


Stable Isotope Analyses 
Our isotope data indicated a broad span of isotopic ori- 


gins for our samples, ranging from xeric to mesic habi- 
tats and different geographic origins (Figure S1). Most 
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individuals had lower wing 67H values than predicted 
considering mean amount-weighted local precipitation 
67H (Figure S1), implying more temperate regions of or- 
igin. In only three cases was the individual wing 87H val- 
ue higher than expected from the local collection site. 
Using an odds ratio of 2:1 or 3:1 resulted in assignment 
of at least 20 out of 29 individuals as immigrants vs. lo- 
cally produced (Table 3). This changed only slightly with 
a more conservative 4:1 odds ratio (19 out of 29). Note 
that these estimates represent minimum numbers of 
immigrants; specimens with 87H values consistent with 
local origin could have originated from a different area 
with a similar local precipitation 67H. 
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Table 3. Stable isotope (57H) results of wing chitin of individual Pantala flavescens taken from subsamples of collections across 
at the sites as indicated. Also indicated are the assignment to immigrant (im) or local (loc) based on given odds ratios as de- 


scribed in the Methods. 











ID Location Latitude Longitude Elevation Wing Precip. Predicted Odds Ratio 

(m) 57H (%o) 67H wing ôH 2:1 3:1 4:1 
(%o) (%o) 
1 India 13.5662 78.499 0 -96.0 -58.7 -58.01 im im im 
2 China 30.4919 119.615 1506 -121.7 -87.0 -108.06 loc loc loc 
3 China 30.4919 119.615 1506 -94.4 -57.0 -108.06 loc loc loc 
4 Japan 35.8896 139.652 12 -145.5 -113.1 -93.50 im im im 
5 Japan 38.5587 140.85 208 -138.4 -105.3 -97.14 im im im 
6 Australia -18.9 145.77 0 -141.9 -109.2 -55.28 im im im 
7 Australia -18.9 145.77 0 -120.1 -85.2 -55.28 im im im 
8 Australia -18.9 145.77 0 -70.6 -30.8 -55.28 loc loc loc 
9 Australia -18.9 145.701 0 -83.2 -44.7 -55.28 im im im 
10 Senegal 14.6433 -12.3363 23 -81.0 -42.3 -67.11 loc loc loc 
11 Senegal 14.6433 -12.3363 23 -100.0 -63.1 -67.11 im im im 
12 Senegal 13.9726 -15.0079 55 -107.5 -71.4 -68.93 im im im 
13 Guinea-Bissau 11.683 -14.7977 14 -92.7 -55.1 -69.84 im im loc 
14 Cruise ship 18.0294 64.1041 (0) -89.7 -51.8 -58.01 im im im 
18 Guyana 6.4833 -58.2167 (0) -126.6 -92.4 -49.82 im im im 
19 Guyana 6.4833 -58.2167 (0) -108.5 -72.5 -49.82 im im im 
20 Hawaii 21.6403 -158.063 5 -76.7 -37.5 -48.91 im im im 
21 Hawaii 21.6403 -158.063 5 -62.8 -22.3 -48.91 loc loc loc 
22 Hawaii 21.6403 -158.063 5 -65.5 -25.2 -48.91 loc loc loc 
23 Hawaii 21.6403 -158.063 5 -77.1 -38 -48.91 im im im 
24 Peru -7.449 -73.9403 295 -88.6 -50.6 -76.21 loc loc loc 
25 Peru -9.8933 -76.32 1981 -96.7 -59.5 -108.06 loc loc loc 
26 Peru -9.8933 -76.32 1981 -89.9 -52.0 -108.06 loc loc loc 
27 Brazil -12.196 -38.9414 250 -137.3 -104.1 -48.00 im im im 
28 Costa Rica 8.6143 -83.5411 304 -115.5 -80.2 -85.31 im im im 
29 USA 35.96 -83.924 260 -130.9 -97.1 -73.48 im im im 
30 USA 30.4642 -84.5953 43 -105.3 -69.0 -58.92 im im im 
31 USA 30.3356 -97.9112 175 -89.2 -51.3 -61.65 im im im 
32 USA 39.8683 -75.1877 4 -116.3 -81.1 -78.94 im im im 
Discussion across five continents for COI and all individuals collect- 


Our genetic data suggest a lack of population structure 
across the global range of P. flavescens, which can be 
attributed to widespread migration. This interpretation 
is supported by our isotope data and by our literature 
review. These findings suggest the existence of a mech- 
anism we call the pantropical Pantala conveyor belt. 
Here, we have compiled the largest P. flavescens hap- 
lotype dataset to date, including for the first time indi- 
viduals from each of the five continents where this spe- 
cies is usually found. Our results demonstrate genetic 
homogeneity across all geographic regions, suggesting 
either a recent shared common ancestry or widespread 
dispersal. With this globally sampled dataset we find 
not only a lack of genetic structuring among regions, 
but also high levels of gene flow, and this suggests 
some level of panmixia. First and foremost, one main 
haplotype is shared among most individuals collected 
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ed across three continents for 16S. Secondly, F, values 
(i.e. estimates of genetic structure within and between 
populations) for COI suggest little population differen- 
tiation and high levels of gene flow between Northern/ 
Southern Hemispheres, Eastern/Western Hemispheres, 
and most continents (Table 2); country dataset results 
were not biologically meaningful, and are likely unre- 
liable given small sample sizes for some countries (cf 
Excoffier, 2007). Thirdly, from our stable isotope analy- 
sis it is clear that a majority of the individuals sampled 
were immigrants to their locations of capture, further 
supporting widespread mixing of P. flavescens. 
Nevertheless, incomplete lineage sorting, or recent 
shared common ancestry, may also be associated with 
such low F., values. In support of this possibility, Tajima’s 
D values (which estimate deviation from neutral expec- 
tations) were negative for all areas for COI, suggesting 
population expansion. These values were statistically 
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significant in the Eastern and Northern hemispheres, 
suggesting less polymorphism than expected under 
neutral mutation and genetic drift, which in turn implies 
that the population is not in equilibrium. This provides 
alternative hypotheses to high gene flow, namely that 
the population may have experienced a genetic bottle- 
neck or may be experiencing purifying selection, either 
of which would also explain the apparent lack of any ge- 
netic structure in P. flavescens across its range. 

K-means clustering (a means of partitioning observa- 
tions into clusters, which in this case partitions individ- 
uals into putative populations) suggests that there may 
be several populations, although none of the “popula- 
tions” identified can be easily interpreted in terms of a 
biological pattern. For example, the clusters contained 
individuals from several countries, continents, and 
hemispheres, across different sampling years, with no 
apparent pattern. Similarly, we did not find any statisti- 
cal support for putative populations defined on the ba- 
sis of various geographic scenarios. F,, values were low- 
est when grouping individuals as Northern/Southern 
Hemisphere or Eastern/Western Hemisphere popula- 
tions. While there is greater apparent genetic structure 
among some continents, this may be biased by unequal 
sampling across regions. 

Migrate-N (a program which models migration us- 
ing genetic datasets) suggests that of two tested sce- 
narios, east-west migration is more probable compared 
to north-south migration; this can be interpreted as 
more common, suggesting the east-west direction is a 
more common migration route than north-south migra- 
tion which makes sense given Pantala flavescens’ use 
of intertropical convergence zone winds. F, values for 
COI recovered between P. flavescens groups were low 
(North-South hemisphere scenario; East-West hemi- 
sphere scenario), which also suggests high gene flow 
among hemispheres (see also Table 2 for additional 
population genetics statistics). The F, values were simi- 
lar to those of the highly migratory Anax junius (micro- 
satellite F = 0.02-0.08, Matthews, 2007), and lower 
than those of somewhat migratory Ischnura hastata 
(Say, 1839) (ND1 F, = 0.105, Lorenzo-Carballa et al., 
2010), or Cordulegaster sarracenia Abbott & Hibbitts, 
2011 (microsatellite F, = 0.423, Abbott et al., 2018). It 
is possible that by using a more rapidly evolving marker, 
such as SNPs or microsatellites, additional structuring 
between areas might be revealed. Alvial et al. (2019), 
however, showed very small differences in Hd, n, and 
number of haplotypes whether COI or microsatellites 
were used for P. flavescens, although with microsatel- 
lites some structure was found between what they char- 
acterised as “continental” and “insular” regions. Their 
continental samples were limited to Central and South 
America. Microsatellites may evolve more rapidly than 
COI (e.g., McDonald & Potts, 1997), and thus if P. flave- 
scens colonized the Americas, most likely from Africa, in 
the relatively recent past, it is possible that such popula- 
tion structure might be apparent in microsatellite data 
(2019) but not in COI data (as demonstrated here). 
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We also note that among the individuals studied 
by Alvial et al. (2019) there was high relatedness only 
among their Rapa Nui samples. This could indicate 
that those samples shared parents, and as a result 
any interpretation of population structure should be 
treated with caution as these are unreliable. Alvial 
et al.’s (2019) samples from Rapa Nui also suggested 
signatures of a recent genetic bottleneck; we do not 
have samples from Rapa Nui, but when their sam- 
ples were included in analyses with our samples we 
found no evidence of genetic structure in COI. If all 
the samples Alvial et al. (2019) collected during one 
trip from Rapa Nui were of closely related individuals, 
that could explain the perhaps artificially high level 
of structure between Rapa Nui and the Americas. Fu- 
ture work should include collection of a time series 
of Pantala from Rapa Nui to determine whether Alvial 
et al.’s (2019) results can be confirmed with samples 
with lower relatedness. Indeed, future work should 
examine more nuclear data to evaluate what biogeo- 
graphical patterns are reflected in the nuclear genome 
of P. flavescens; there can be genetic discordance be- 
tween gene fragments due to phenomena such as in- 
complete lineage sorting and introgression, while pat- 
terns of events such as hybridization and isolation can 
lead to biogeographical discordance (Toews & Brels- 
ford, 2012). In addition to augmenting knowledge of 
Rapa Nui, future efforts should also aim to increase 
genetic sampling across the continent of Africa. More 
generally, we encourage further study with more and 
larger samples, and a range of genetic markers. We 
note in particular that for COI (for which at least 71 
haplotypes have been identified) very much larger 
sample sizes than used here would be desirable, and 
also that the study of faster-evolving genetic markers 
than COI may be needed to fully elucidate the popula- 
tion structure of P. flavescens. 

Our stable isotope data support the general percep- 
tion that P. flavescens is a highly vagile species. The 
majority of individuals investigated here were classed 
as immigrants, which means we collected individuals 
from areas other than their natal regions regardless 
of where or when we sampled. Stable isotope analy- 
sis has previously shown that P. flavescens collected in 
Maldives were all immigrants (Hobson et al., 2012a), 
while those collected in Japan were mostly immigrants 
(Hobson et al., 2021). The fact that we observed only 
small numbers of potentially local individuals in our 
samples confirms that most adults disperse soon after 
emergence (e.g. Corbet, 1964; Kumar, 1972). Increased 
isotope sampling in areas like the Andes, where we had 
only potentially local individuals, and particularly on 
isolated islands (coupled with paternity studies) could 
assist with establishing whether there are exceptions to 
the patterns suggested by our results here, (i.e. wheth- 
er there are any resident non-migratory populations 
in other parts of the range of P. flavescens). However, 
even in Hawaii, two out of four sampled individuals 
were classed as immigrants. 
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Finally, a literature review confirmed large-scale 
movements of P. flavescens between locations where it 
is commonly found (text $1, Table $3; Figure S2). There 
are several possible pathways of population exchange, 
often as large swarms, between continents and even 
across large stretches of open ocean (text S1). Most of 
these observations were made in association with par- 
ticular wind patterns such as the Inter-Tropical Conver- 
gence Zone (ITCZ) rain fronts, but also associated with 
typhoons and cyclones. Noted dispersal events imply 
that individuals are able to cover distances of several 
thousand kilometers even over oceans and seas, e.g., 
Anderson, 2009; Drake et al., 1995; Sparrow et al., 
2020). This supports the idea that migration and dis- 
persal may be wind-driven or wind-supported. We as- 
sume that seasonal appearances at off-shore islands can 
provide clues for reconstructing cross-ocean migration 
routes (cf. Anderson, 2009). One such migration route 
has been described by Anderson (2009) crossing the In- 
dian Ocean between India and Africa via the Maldives 
and Seychelles Islands, and there is additional evidence 
for cross-Indian Ocean travels (text S1); Hobson et al. 
(2021) also describe oceanic crossing by individuals col- 
lected on Japanese islands. As many Pacific islands are 
also visited by P. flavescens, there is support for the 
idea of intercontinental pathways via the Pacific Ocean. 
Whether P. flavescens can regularly cross the eastern Pa- 
cific barrier remains unclear. The most likely route con- 
necting the Eastern and Western Hemispheres seems 
to be via the central Atlantic Ocean from West Africa, 
where individuals would be supported by trade winds 
blowing towards northern South America and the Carib- 
bean. This dispersal path has been demonstrated for two 
other migratory Odonata of African origin, Anax ephip- 
piger (Burmeister, 1839) and Tramea basilaris (Palisot 
de Beauvois, 1807) (Lambret et al., 2013; Hedlund et 
al., 2020; Meurgey & Poiron, 2012) as well as the desert 
locust Schistocerca gregaria (Forsskal, 1775) (Lovejoy et 
al., 2006). Seasonal winds may also assist the movement 
of swarms of P. flavescens as far away from the tropics as 
about 52°N in East Asia (Borisov, 2012; Cao et al., 2015), 
North America (May, 2013) and, recently, also in parts of 
Europe (Buczyński et al., 2014); in these high latitude re- 
gions they probably do not survive the winter (Borisov, 
2012; Borisov & Malikova, 2019; May, 2013). In Asia and 
Africa, movements seem to be mainly associated with 
the seasonal shifts of the ITCZ and the associated rain- 
bearing winds, while reasons for massive overshoot of 
the ITCZ in North America are not yet fully understood 
(May, 2013). In the Southern Hemisphere most larger 
landmasses except for Antarctica and southern Patago- 
nia regularly receive migratory swarms of Pantala (al- 
though P. flavescens has been found at least as far south 
as 37°S in Argentina; del Palacio et al., 2017). Our con- 
clusion, based on observations compiled from literature 
(supplementary material), is that the extraordinary dis- 
persal abilities of Pantala flavescens coupled with sea- 
sonal wind patterns, support the observed global pat- 
tern of gene flow and migration. 
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Physiological and morphological features that may 
allow individuals to efficiently migrate long distances 
likely involve their ability to efficiently use lipids as fuel 
rather than carbohydrates (e.g., Kallapur & George, 
1973; it has been suggested that dragonfly fat reserves 
can sustain eight hours of flapping flight), as well as dis- 
tinct wing morphology (e.g., Alvial et al., 2019; Moore, 
1993; Sacchi & Hardersen, 2013; Sudrez-Tovar & Sar- 
miento, 2016; Outomuro & Johansson, 2019; Zhao et 
al., 2012), and unique thoracic musculature (Bäumler 
et al., 2018). Further, a major precondition (or trait) for 
being an obligate dragonfly migrant may be the abil- 
ity to develop successfully in ephemeral freshwater 
habitats, which are usually formed after seasonal rains 
(Corbet, 1999; Suhling et al., 2019). This seasonality in 
breeding habitat is directly linked to the propensity of 
P. flavescens to follow the ITCZ. 


Conclusions 


Using COI data from samples taken across most of P. fla- 
vescens’ range, we found high gene flow and a lack of 
genetic structure, with one main haplotype shared 
globally. This result is supported by isotope data which 
suggest frequent long-distance migration. We also com- 
piled extensive observational evidence from the litera- 
ture which supports these findings. We conclude that 
there is full genetic exchange among all continents, i.e., 
a globally panmictic population of P. flavescens seems 
possible, even if very isolated islands may be excluded 
from regular migration. The mechanism supporting pan- 
mixia is regular long-range migration within and between 
continents, even across oceans. Trans-oceanic migration 
sustained by trade winds and additional storm events 
explain the ability of individuals to reach distant islands 
and to travel between continents, while active flight and 
movement of seasonal rains in relation to the ITCZ can 
explain patterns observed over continents. All evidence 
points at a mechanism we call the pantropical Pantala 
conveyor belt, a mechanism which makes it possible for 
an insect with excellent flying ability and rapid develop- 
ment to achieve some level of global panmixia. 
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Supplementary material 


Supplementary Figure S1. Patterns of wing isotopes (a) world map 
generated from Google Earth (Map data ©2020 Google, earth. 
google.com/web/, United States) with wing isotope ranges in- 
dicated for subsample of 32 specimens, map generated with 
Google Earth (b) scatterplot of wing isotope values versus pre- 
dicted wing isotope values, the dashed line depicts the 1:1 situa- 
tion, the shaded area the confidence interval. 


Supplementary Figure S2. Map from Google Maps (Map data 
©2020 Google, earth.google.com/web/, United States) gener- 
ated in Excel, compiling data from Tables 3 and S3. Intertropical 
convergence zone (ITCZ) displayed as a line. TOM = transoceanic 
migration, CMS = coastal migratory swarms, TMS = terrestrial mi- 
grating swarms, LDD = long distance dispersal events to formerly 
uninhabited areas. 


Supplementary Figure S3. PopART minimum spanning network; Da- 
tasets: Cytochrome B (N = 77), Histone 3 (N = 4), 16S (N = 16), 
and 18S (N = 4). Haplotype networks coloured based on sam- 
pling country unless only one country of origin for the samples, 
in which case the networks are coloured black. 
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